home *** CD-ROM | disk | FTP | other *** search
/ Resource for Source: C/C++ / Resource for Source - C-C++.iso / codelib7 / v_09_09 / 9n09016a < prev    next >
Encoding:
Text File  |  1995-11-01  |  1.6 KB  |  59 lines
  1. /* _Exp function */
  2. #include "xmath.h"
  3.  
  4. /* coefficients, after Cody & Waite, Chapter 6 */
  5. static const double p[3] = {
  6.     0.31555192765684646356e-4,
  7.     0.75753180159422776666e-2,
  8.     0.25000000000000000000e+0};
  9. static const double q[4] = {
  10.     0.75104028399870046114e-6,
  11.     0.63121894374398503557e-3,
  12.     0.56817302698551221787e-1,
  13.     0.50000000000000000000e+0};
  14. static const double c1 = {22713.0 / 32768.0};
  15. static const double c2 = {1.428606820309417232e-6};
  16. static const double hugexp = {(double)HUGE_EXP};
  17. static const double invln2 =
  18.     {1.4426950408889634074};
  19.  
  20. short _Exp(double *px, short eoff)
  21.     {   /* compute e^(*px)*2^eoff, x finite */
  22.     int neg;
  23.  
  24.     if (*px < 0)
  25.         *px = -*px, neg = 1;
  26.     else
  27.         neg = 0;
  28.     if (hugexp < *px)
  29.         {   /* certain underflow or overflow */
  30.         *px = neg ? 0.0 : _Inf._D;
  31.         return (neg ? 0 : INF);
  32.         }
  33.     else
  34.         {   /* xexp won't overflow */
  35.         double g = *px * invln2;
  36.         short xexp = (short)(g + 0.5);
  37.  
  38.         g = (double)xexp;
  39.         g = (*px - g * c1) - g * c2;
  40.         if (-_Rteps._D < g && g < _Rteps._D)
  41.             *px = 1.0;
  42.         else
  43.             {   /* g*g worth computing */
  44.             const double y = g * g;
  45.  
  46.             g *= (p[0] * y + p[1]) * y + p[2];
  47.             *px = 0.5 + g
  48.                 / (((q[0] * y + q[1]) * y + q[2]) * y
  49.                 + q[3] - g);
  50.             ++xexp;
  51.             }
  52.         if (neg)
  53.             *px = 1.0 / *px, xexp = -xexp;
  54.         return (_Dscale(px, eoff + xexp));
  55.         }
  56.     }
  57. /*End of File */
  58.  
  59.